Mon. Not. R. Astron. Soc. 000, 000-000 (0000) Printed 3 February 2008 (MN lAT^ style file v2.2) 



On time-dependent orbital complexity in gravitational 
N-body simulations 

N. T. C. M. Boily^ and S. Portegies Zwart^'^ 

^ Observatoire Astronomique, Universite de Strasbourg and CNRS UMR 7550, 11 rue de I'Universite, 67000 Strasbourg, France 
^ Sterrenkundig Instituut 'Anton Pannekoek' , University of Amsterdam, Kruislaan 403, 1098 SJ Amsterdam, The Netherlands 
^Section Computational Science, University of Amsterdam, Kruislaan 4-03, 1098 SJ Amsterdam, The Netherlands 



oo 

o 

O 
(N 

D 
(N 



6 



> 

cn 
cn 

O 

00 

O 

> 

X 



3 February 2008 



ABSTRACT 

We implement an efficient method to quantify time-dependent orbital complexity in 
gravitational iV-body simulations. The technique, which we name DWaTIM, is based 
on a discrete wavelet transform of velocity orbital time series. The wavelet power- 
spectrum is used to measure trends in complexity continuously in time. We apply the 
method to the test cases iV = 3 Pythagorean- and a perturbed N = 5 Caledonian 
configurations. The method recovers the well-known time-dependent complexity of 
the dynamics in these small-iV problems. We then apply the technique to an equal- 
mass coUisional N = 256 body simulation ran through core-collapse. We find that a 
majority of stars evolve on relatively complex orbits up to the time when the first 
hard binary forms, whereas after core-collapse, less complex orbits are found on the 
whole as a result of expanding mass shells. 
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1 INTRODUCTION 

Computer simulations of the gravitational A^'-body problem 
aim to solve the set of SA'^ second order ordinary differential 
equations 



E 



mi{ri 



\rr - ri\ 



,N). (1) 



Here G is the gravitational constant, Fi and f'i denote the 
force and the acceleration exerted by the A'^ — 1 particles 
I on particle i of mass rrii at position ri. A system with 
A' > 2 is in general chaotic, i.e., the solution to Eq. ([TJ is 
known to be sensitive to the initial conditions. Miller (1964) 
was the first to show that individual orbits calculated from 
neighbouring initial configurations diverge on a short time- 
scale proportional to 4^/A', where (, is the mean time be- 
tween two subsequent close encounters of a particle. Even 
when computed with double precision arithmetic, the dif- 
ferences between the two integrations become comparable 
to the characteristic length- and velocity-scales of the clus- 
ter within only a few crossing times. Thus one may argue 
that the positions and velocities obtained from an A'^-body 
simulation are a fair rendition of the problem at hand in 



Tremaine 19921. Goodman, Heggie & Hut (19931 identify 



three mechanisms responsible for the exponential growth of 
the distance between nearby trajectories. Aside from numer- 
ical error^ they point out that the exponential instability of 
the solutions also arises through a fluctuating mean gravita- 
tional field, or through inherently chaotic orbits in a static 
field. However, for systems in dynamical equilibrium and 
with near-spherical symmetry (such as e.g. globular clus- 
ters), they conclude that the principal mechanism respon- 
sible for chaos is the cumulative effect of near neighbour 
interactions, i.e. two-body encounters. 

The exponential divergence of the A''-body problem of 
Eq. ([T]) has been studied by several authors for sys tems with 
up to A'^ = 10^ particles ( Kandrup fc Smith|l991 Kandrup 



& Sideris 2001 Hemsendorf & Merritt 20021. The debate 



whether Miller's instability is formally caused by chaos is 
still on-going (see e.g., [Kandrup fc Sideris|[2003} [Helmi fc| 
G6mez|2007| . The assumed chaoticity of such systems may 



be evaluated by computing the variational equations associ- 
ated to Eq. ([T| ( |Miller|1971| ) and by retrieving indicators of 
chaos such as, for instance, Lyapunov characteristic numbers 
(Bcnettin, Galgani fc Strelcyn 19761. Other indicators of 



a statistical sense only ( [Aarseth fc Lecar]|l975| |Quinlan fc 



chaos commonly fou nd in the literature are e.g., the relativ e 
Lyapunov indicator ( Sandor, Erdi fc Efthymiopoulos|2000 1, 
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^ A round-off error at one time step can be considered as a change 
in the initial conditions for the next time step. 
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the mean exponential divergence of nearby orbits ( Cincotta 
& Simo 2000 1 or the small- alignment index ( Skokos||2001 1 
All these indicators accurately quantify the exponential di 
vergence of nearby trajectories. 

In this work, we propose to investigate time-dependent 
orbital complexity in an A'^-body simulation. We define or- 
bital complexity to be a measure of the richness and non- 
triviality of the frequency spectrum of an orbit at a given 
time t. The possible connection between orbital chaos and 
orbital complexity has been pointed out by |Kandrup, EcF] 
stein & Bradley (19971. In that work, the authors find a 



strong correlation between their measure of orbital com- 
plexity and short time Lyapunov exponents. The notion of 
complexity provides information about the orbital content 
of a gravitational system governed by Eq. ([T]) that is com- 
plementary to the classical indicators of chaos mentioned 
above. For instance, the concept of complexity is exploited 



by Sideris & Kandrup ( 2002 1 to study the continuum limit 



in the case of large-A'^ simulations. It is used to compare 
the orbital behavior between particles evolving in smooth 
potentials and bodies orbiting in the corresponding frozen 
Ai'-body configurations. The notion of complexity here allows 
to contrast the discreteness effects of the A'^-body configu- 
ration to the orbital evolution obtained in the smooth case. 
However, whereas a global measure of orbital complexity 
has been implemented in several works (see e.g., Kandrup 
et al. 1997), a method to measure the impact of instanta- 
neous changes in orbital complexity has not. A trajectory 
computed from Eq. ([TJ can show multiple, qualitatively dis- 
tinct regimes in time. For instance, an orbit may display 
arcs of relatively smooth motion such as e.g., an unper- 
turbed parabolic or hyperbolic orbit. At other times, the 
body may be gravitationally bound in a binary system with 
a particle of approximately the same mass. Likewise, the 
body may be temporarily trapped in a complicated higher- 
order resonance, orbiting about a massive central body. All 
these states of motion can be identified in time by a suitable 
measure of complexity. The goal of this paper is precisely to 
discuss such time-resolved complexity by introducing a ded- 
icated tool for complexity evaluation. The issue of formally 
relating complexity to chaos will not be addressed here. 

Classical spectral methods and Fourier transform based 



techniques (see Laskar 1993 Carpintero & Aguilar 1998 
Valluri fc Merritt|1998[ ) are of short execution time and able 



to provide an accurate frequency domain representation of a 
given orbit. A global measure of complexity can be retrieved 
by such an approach. However, the Fourier transform suffers 
from the problem of losing any time-dependent information 
on the motion of the particle. In this paper, we present an 
accurate, easy-to-implement and time-resolved complexity- 
detection tool for individual orbits in A'^-body simulations. A 
state of dynamical equilibrium is assumed throughout. We 
implement a specially-adapted method of time series analy- 
sis that is based on a discrete wavelet transform. The method 
post-processes the orbital data of a simulation and provides 
a time-resolved measure of complexity. The present work 
mainly focuses on a detailed description of the method and 
is structured as follows. In Section 2 we introduce the identi- 
fication technique for time-dependent complexity. Section 3 
presents applications to A'^ = 2, 3 and 5-body problems. The 
analysis of a low-resolution A'^ — 256 equal-mass Plummer 



model is presented in Section 4 and a brief discussion of the 
results is given in Section 5. 



2 METHOD, TESTS AND VALIDATION 

This section presents and tests a technique to quantify time- 
dependent orbital complexity. In §2.1| we present the dis- 
crete wavelet transform (hereafter DWaT) as an efficient 
tool to determine the base frequencies of an orbit in a time- 
resolved manner. The indicator of complexity obtained from 
the DWaT, namely the discrete wavelet transform informa- 
tion measure (DWaTIM), is defined in jp^ 



2.1 Wavelet transforms 

2.1.1 Definitions 

The complexity of motion is analyzed by means of a wavelet 



transform (see also Gemmeke, Portegies Zwart & Kruip 
|2008[ ). The use of wavelet transforms has found applications 
in a wide range of domains such as seismology, financial 
time series processing or medical electrocardiogram studies; 
see |Hubbard| p998| ) for further references. A wavelet trans- 
form provides a time-frequency representation of a time se- 
ries f{tq) by fitting a wavelet 5' to a set of points {tq}. 
Whereas Fourier-based methods decompose a signal into in- 
finite sine and cosine functions, effectively losing information 
at individual times tq , the wavelet transform offers precise lo- 
calization in both the frequency- and time-domain. Wavelet 
transforms remain band-limited however; they are made up 
of not one but a limited range of several frequencies. 

A wavelet family ^'a,T is defined by the set of elemen- 
tal functions generated by scaling and translating a mother 
wavelet "^{t): 



*„,,(*) = |a| 2* 



(2) 



where a represents the scale variable and r the translation 
variable (a, r £ K, a 5^ 0). 

The continuous wavelet transform (CWT) is defined as 
the correlation between a signal S(t) G L^(R) (the space of 
square summable functions) and the wavelet family ^'a,T(i) 



for each a and r (see e.g., Daubechies||1992 1 



{Sit),^a.r) = \a\ 



S{t)-^ 



dt. 



(3) 



Here {S, *I'a,r) denotes the wavelet coefficients and ^ is the 
complex conjugate of 'if. Equation ([3| can be inverted to 
reconstruct the original time series. The CWT is known to 
produce a large amount of wavelet coefficients which im- 
plies considerable CPU execution times (see e.g 



Samar et 



al. 19991. In addition, the information the CWT displays 



at closely spaced scales or at closely spaced time points is 
highly correlated and thus unnecessarily redundant. 

For these reasons we instead compute a discrete wavelet 
transform (DWaT). The DWaT offers a highly efiicient 
wavelet representation that can be implemented with a sim- 
ple recursive filter scheme ( Mallat|[l999 Daubechies|| 1 992 1 . 
Unlike the numerical CWT implementation which easily 
produces more than 10^ coefficients for a single orbital time 
series of Q = 8192 data points, the DWaT only produces 
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as many coefficients as there are samples in the time series, 
i.e. Q. This property of the DWaT of avoiding redundant 
wavelet coefficients serves in defining a proper measure of 
complexity, as we will show in §2.2| For a given choice of 
the mother wavelet function '5(4) and for the discrete set of 
parameters aj = 2-' and rj,fc — 2-'k (j, k £ Z), the wavelet 
family 



t - T. 



3,k 



= 2"3*(2"^f-fc) 



(4) 



defines an orthonormal basis of L^{W). The time series f{tq) 
is sampled at Q = 2"^ ( J £ N) constant time intervals of size 

A = tq + l — tq 

(see e.g.. Rosso et al.|2006 Equation [36]) 



(5) 



J 

For simplicity we set A = 1 throughout ^ Here f{t) is the 
reconstructed signal and J = log2 Q is the number of scales 
over which the time series f{tq) is analyzed. The DWaT 
coefficients (/(t,), ^'j,fc) discrete can be understood as a rep- 
resentation of the wavelet power spectrum (or, energy) at 
scale j and time tq, associated to the time series /(ig). They 
represent the local residual errors between successive signal 
approximations at scales j and j — 1. In what follows, we set 



(6) 



for a more convenient notation of the DWaT coefficients. 
As mentioned earlier, there are Q coefficients Pj{k) and the 
number of coefficients computed for resolution level j is 2-' "^ . 
The frequency band over which the Pjik) are computed is 
limited by the frequency 1/(QA) (scale j = 1) in the low- 
frequency domain and by the Shannon-Nyquist critical fre- 
quency /cr = 1/(2 A) in the high-frequency domain (scale 
j = J; see also|Press et al.|2002[ §12.1 and §13.10). We refer 
to |Samar et al.| ( |1999[ ) for further details about the wavelet 
representation. 

We use bi-orthogonal cubic spline functions as mother 
wavelets ( [Cohen, Daubechies fc Feauveau| |1992| case 
{N,N} = {3,9} in their Table 6.1), 



*(t) = ^»,e3(2t-/), 



(7) 



where the (jj's are known as basic spline coefficients and 



-(f- 1/2)2 + 3/4, 

(t-2)V2, 

0, 



-1 < i < 

< t < 1 

1 < i < 2 
otherwise. 



(8) 



This choice is motivated by three arguments. First and 
most importantly, spline functions provide an excellent 
time-frequency localization when compared to other mother 
wavelet candidates ( Ahuja, Lertrattanapanich fc Bose|2005 



Unser|1999[ ). Instantaneous changes in the dynamics are ac- 
curately singled out by the DWaT. We further stress this 
point in §2.3| Second, the use of splines is computationally 
inexpensive ( [Thevenaz, Blue fc Unser||2000[ ) and provides 
further desirable properties such as e.g., compact support 
and smoothness. (For an exhaustive discussion on spline in- 
terpolation, see |Unser|1999l [Thevenaz, Blue fc Unser|2000[ ) 
Finally, the use of a bi-orthogonal spline mother wavelet also 
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Figure 1. Border effects of the discrete wavelet transform 
(DWaT) of the two sinusoids y\ and y2 defined in Eq. |10| From 
top to bottom: (a) time series of y\ and y2 (the dashed line is 
y2), (b) DWaT of yi, (c) DWaT of y2 with full border effects, (d) 
DWaT of y2 with reduced border effects. 



implies reduced border effects, an undesired artifact of the 
wavelet transform algorithm (see ^2.1.2| below). 
We analyze the velocity time series 



f(tq) = Vai(tq) 



(a = x,y,z), 



(9) 



i.e. the v^i , Vy- and v^^ velocity components of each particle 
i {i = 1, ...,N). We do not use the information available on 
the positions of the bodies since there may be important 
differences in magnitude between the beginning and the end 
of these time series. Positional information is then likely to 



produce pronounced DWaT border effects (see S 2.1.2 1 



2.1.2 Border effects 

Border effects are an artifact of the wavelet transform algo- 
rithm which enforces cyclical boundary conditions on the 
data vector ( [Lo Presti fc 01mo|[T996| [Press et aL][2002| 
§13.10). Such border effects depend on the values held by 
the two end points of the data vector, and may become im- 
portant when both ends of the data set differ greatly. In the 
following, we aim to quantify the extent to which the diag- 
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nostics become erroneous due to these edge effects. 
Figure [T][a) shows the two sinusoids 



yi 



2/2 



32 Stt 

cos ivr X q H 

8192^ 2 



/ o 15/2 
cos I Itt X q 



8192 



(10) 



sampled at Q = 8192 intervals (g — 1,...,Q). The inte- 
gration time is commensurate to the periodicity of signal 
2/1. This is not the case for signal 1/2. Figure [TJb) and (c) 
show their respective DWaT scalograms. The scalogram is a 
grey-shaded representation of the DWaT coefficients Pj(k). 
The darkest shade is for the largest values |Pj(fc)|; white 
means Pj{k) — 0. For both sinusoids a maximum intensity 
is obtained for the scale that represents the base period Ilj 
of the respective signal. These are scales j = 6 at period 
= 2^-^ ■ QA = 2"^ ■ 8192 = 256 for yi and j = i for 
2/2 . The region in which the discrete wavelet transform of j/2 
suffers from border effects is known as the cone of influence 
i Moortel, Munday fc Hood|[2004l l . The cone of influence is 
clearly visible at both edges on Fig. [ijc). On the left-hand 
edge of Fig.[TJc), for example, the DWaT gives an artifl- 
cially high excited mode at scale j — 5 up to t = 1024. The 
same artifact is also found at scales of higher frequencies, 
although the magnitude of the effect then diminishes and is 
barely visible on the left-hand side of Fig.[TJc). In addition, 
scale J = 4 is wrongly excited in the intervals < t < 3072 
and 7168 < t < 8192. The size of the cone of influence de- 
pends on the choice of the mother wavelet and is especially 
significant when the spectrum of the signal contains low fre- 
quencies, when the period compares to the overall duration 
of the time series (Moortel et al. 2004). 

Figure [T][d) shows a situation where the border effects 
have been reduced by extending the integration time to twice 
the original intervaF] Doubling the time interval of analy- 
sis implies an increased ratio between the signal length and 
its periodicity. This allows for a higher resolution in the fre- 
quency domain (Carmona et al. 1998). Furthermore, a larger 
portion of the DWaT remains unaffected by the border dis- 
continuity and thus a larger amount of reliable information 
can be retrieved. In particular, the left-hand and the right- 
hand border effects are diminished on Fig.[l][d): for instance, 
the wrongly excited mode at j = 5 between < t < 1024 
in Fig.[T][c) has been reduced by « 82% and the scale j = 4 
artifacts of the right-hand side of that figure have vanished 
completely. In the remainder of this paper, we apply the 
technique of time series extension whenever the edge effects 
appear significant. To take into account the remaining ar- 
tifacts at the beginning of the time series, we furthermore 
analyze the first ~ 2048 data points (equivalent to the first 
~ 25% of the signal) of the DWaT analysis with particu- 
lar caution. Although Fig.[T][d) indicates that border effects 
may influence more than the first 25% of the 2/2 analysis, 
we found that this is in general not the case for signals pro- 
duced by A'^-body orbits (see j j2.2.2| and !j3|. In what follows, 
we refer to these 25% of the analysis as potentially biased 
due to border effects. 



^ Alternately, the j/2 data vector may also be padded with zeros 
until twice the original time interval is reached. 



2.2 Measures of complexity 

Our goal is to obtain from the DWaT a quantitative estimate 
of the time-dependent complexity of the frequency spectrum 
of a trajectory. Exploiting some notions of information the- 
ory, we here present the discrete wavelet transform informa- 
tion measure (DWaTIM) as an efflcient indicator for com- 
plexity. In what follows, we provide a succinct overview of 
the concept. We follow closely the approach of [Rosso et al.| 
(2006) and [Martin, Plastino fc Ro sso' f2006|. We refer the 
reader to these works for a more extended discussion. 



2.2.1 Discrete wavelet transform information measure 
(DWaTIM) 

For a chosen time window of size kA (where k is an arbitrary 
integer) we compute the wavelet energy at each resolution 
level j, 



k = l 

and the total wavelet energy. 



(11) 



(12) 



to obtain the so-called relative wavelet energy pj = Ej/E at 
scale j. The DWaT then provides a probability distribution 



V^{P,}, (j = l,...,J), 



(13) 



which weighs the base frequency j in the reconstruc- 
tion of the original signal (tg). By definition we have 

j:Up^= 1- 

The amount of disorder present at time t in an orbital 
time series can be quantified by determining the information 



needed to describe the orbit at that time (Shannon 19481 



An information measure (hereafter IM) can be seen as a 
quantity that describes the characteristics of the time-scale 



probability distribution V of Eq. (13l (Rosso et al. 20061 



The IM gives the amount of information required per time 
unit to specify the state of the system up to a given accu- 
racy. The IM we use in this work is the Shannon entropy 
( Shannon|[1948 Quian Quiroga, Rosso fc Basar[[1999 Sello 
2003f 



(14) 



where = 1 is an arbitrary numerical constant. A minimum 
of information entropy min (Vys['P]) = is obtained if pj — 1 
for some scale j and for all the remaining scales. This 
situation only occurs for the ordered dynamics of a periodic 
orbit with a single base frequency. Likewise, the state of 
highest complexity is obtained for the case of a white noise 
signal. For such an orbit, the entire band of base frequencies 
is sampled by the DWaT in equal proportions, and V is 
characterized by the uniform probability distribution 



Vu = 



(15) 



In this case the IM is max(M/s[P]) = Ws[Vu] = log2 J. For 
a given velocity component u^. (a — x,y,z; i = 1, A^) we 
define the DWaTIM to be the measure 
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Figure 2. Complexity measure for the sinusoids yi and j/2 of 
Fig.^ The thick solid line confounded with the x-axis shows the 
DWaTIM Hy^ as obtained for yi from Fig.[l];b). The upper solid 
line shows the DWaTIM Hy^ obtained for y2 from Fig.[Tj[d) in 
which border effects where reduced by time series extension. The 



border effects have not been treated (see Fig.[TJc]). 



Ha- 



( Ws[V] 



(16) 



i.e. the Shannon entropy for velocity component a of particle 
i normalized to the interval [0, 1]. The generalized DWaTIM 
of particle i is then obtained by taking the arithmetic mean 
over the 3 components, 

+ Hy. + Hz- 



(17) 



Finally, the overall DWaTIM Ttot of the system is computed 
by averaging over all particles i, 



Ttot 



(18) 



iv ■ 

In what follows we will argue that an increase in DWaTIM 
correlates with an increase in the complexity of the under- 
lying orbital dynamics at that instant. 

2.2.2 Complexity of a sinusoid 

The temporal evolution of the DWaTIM indicator is com- 
puted by subdividing the input signal in non-overlapping 
time windows of size kA. We use k — 2. Figure |2] shows 
our complexity analysis for the sinusoids yi and y2 of §2.1.2| 
The DWaTIM of yi, Hy^ is shown by the lower thick solid 
line. As expected for a stationary signal with a unique base 
frequency, Hy-^ ~ 6.13 • 10~*; the line almost coincides with 
the J/ = axis. In what follows, this value can be considered 
as indicating zero complexity. 

We now discuss the pathological case of signal y2 al- 
ready mentioned in §2.1.2| The accuracy of the DWaTIM 
depends on the frequency of the sinusoid. The DWaT has 



a better resolution at high frequencies (Samar et al. 19991 



Conversely, it is more difficult to obtain reliable measures of 
complexity for signals of low frequencies. In addition, bor- 
der effects are the most important for such time series (see 

^ For a more convenient notation, we denote the DWaTIMs ob- 
tained for the test cases yp (P = 1, 2, 7, //) presented in section 
^hyHy^. 



\ 2.1.2 I. The time series y2{tq) is a long-period, non-cyclical, 
unique-frequency signal. The DWaTIM Hy^ is shown by the 
upper solid line on Fig.[2] The line depicts the DWaTIM ob- 
tained from the DWaT of Fig.[l][d), i.e. for the case where 
border effects have been reduced by time series extension. 
(Note that the border effects for that case still influence the 
result at all times.) For the first 25% of the diagnostic the ef- 
fects are dramatic: we find complexity values of up to ~ 0.4, 
where « is expected. The curve Hy2 approaches zero for 
t > 6500 only. Our method of border effect reduction thus 
brings little improvement to the quality of the diagnostic for 
that signal (compared with no corrections at all, cf. dotted 
line on Fig.[2] obtained from Fig.[ljc]). However, one should 
bear in mind that the DWaT approach stems from the idea 
of performing a mtt/ti-frequency analysis of signals that are 
well-resolved over the time interval of interest. It is difficult 
to analyze with the DWaT non-cyclic single-frequency sig- 



nals of long period such as e.g., signal yu (see also Moortel, 
Munday fc Hood||2004| |Samar et ar||1999[ ). Our analysis of 
such sinusoids exposed the limits of the DWaTIM method; 
we here obtained complexity results that were widely bi- 
ased by border effects. However, in this work we aim to 
discuss the more generally encountered multi-frequency sig- 
nals describing orbital motion in the gravitational A^-body 
problem. For this type of signal, we found that the complex- 
ity analysis yields in general consistent results whenever the 
lowest frequency component of the time series is resolved 
over at least ~ 8 complete oscillations. Border effects then 
play a minor role. (For a further discussion of this, see the 
results obtained for the multi-frequency signals of §2.3[ §3.1| 



or i3.2 



2.3 Toy Models 

We illustrate the capability of the method to capture time- 
dependent changes in the dynamics with two toy models. 
Let us consider the two time series yi{tq) and yii(tq) {q — 
1, 8192; tq — qA = q) constructed as follows: 



and 



yii{tq) = < 



VA, 



VA, 

yc, 

Vb, 



<q< 4096, 
4096 < g < 8192, 



< g < 1024, 
1024 <q< 3072, 
3072 <q< 4096, 
4096 <q< 6144, 
6144 <q< 8192. 



(19) 



(20) 



The series yi and yu of Eqs. ( |19[ ) and ( |20[ ) are made up 
of 2 and 5 different dynamical regimes, respectively. Each 
is a different linear combination of Fourier components and 
white noise: 



VA 



Vb = 



sin 27r X 



sin 27r x 



-I- sin 27r x 



32 
8192 

32 
8192 

128 



q + TV 

q + -K 
57r 



8192'"^^ 4 



(21) 



(22) 



yc = 



sin 27r x 



32 
8192' 



g M- sin 27r X 



128 tt' 
8192*^^ 1 , 
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Figure 3. Analysis of the time-dependent complexity of the toy 
model defined by Eq. |19[ From top to bottom: (a) time series, (b) 
discrete wavelet transform, (c) discrete wavelet transform infor- 
mation measure Hy , . 




Figure 4. Analysis of the time-dependent complexity of the toy 
model defined by Eq. |20[ From top to bottom: (a) time series, (b) 
discrete wavelet transform, (c) discrete wavelet transform infor- 
mation measure Hy , , . 
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VE 



yc + Gsiq), 
Si (9). 



(23) 

(24) 
(25) 



Here Qaiq) and Gi{q) are Gaussian white noise functions 
of amplitude unity, zero mean and an absolute variance of 
3 and 1, respectively. The complexity increases from one 
regime to the next: as we go from regime tja, through regime 
ys up to yc, we double the number of frequency components 
each time. We then construct i/d by adding the noise func- 
tion Gsiq) to the yc signal. Finally, regime j/_e consists of 
random numbers only and defines our signal of maximum 
complexity, i.e. with no underlying periodic signal. In what 
follows we study the performance of the DWaTlM to cap- 
ture the time-dependent complexity of the toy models yi 
and yii. 

The complexity analysis for toy model yi of Eq. ( |19[ ) 
is represented in Fig.|3] Figure |3]^a) shows the time series 
yi{tq). Until t — 4096, the signal has the same unique base 
frequency as the sinusoid yi of §2.1.2| However, at f = 4096 
the behavior of the signal changes abruptly. A second com- 
ponent, with a frequency 4 times higher than the primary 
component, is added to the signal. Figures [3|b) and (c) show 
the associated complexity analysis. The DWaT scalogram of 
Fig.jsjb) illustrates the power spectrum of the base frequen- 
cies 1 /Ilj as a function of time. The variation of complexity 
with time of signal yi is clearly identifiable, as the increase 



in complexity a,t t — 4096 is well recovered. The modes rep- 
resented by scales j — 8 and j = 9 show a response between 
4096 < t < 8192. The related DWaTIM complexity measure 
Hyj is plotted in Fig.[3];c). Between 4096 < t < 4200 the 
complexity rises from Hy^ ~ 6 ■ 10"* to an average value of 



Hy 



0.16 for t > 4200. We note that in the time interval 



4096 < t < 4200, Hy-^ rises progressively. The size of this 
transient interval corresponds to about one half of the pe- 
riod of yA which is 256 units. This is the longest of the two 
periods of the Fourier components involved in the transition 
yA — > yB- This sets a limit on the time-localization of the 
transition. In §3.2| we will discuss in detail the latency of 
transitions in the dynamics measured by the DWaTIM in 
an application to small-A'' problems. We also remark that 
border effects between < t < 200 affect the measure for 
less than 3% of the signal. 

Figure J4] shows the results for the second toy model 
yii of Eq. ( |20[ ). The related time series is shown in Fig.|4|a) 
and the results for the complexity analysis are represented 
in Fig.[4]^b) and (c). The time-dependent subtleties of the 
frequency spectrum are depicted by the DWaT in Fig.[4|b). 
For instance, the drop in complexity at t ~ 3072, where we 
remove the two incommensurate base frequencies of regime 
yc and switch to regime j/s, is recovered. The modes j — 10 
and 11, excited between 1024 < t < 3072 (the j = 11 re- 
sponse during that interval is hardly visible on Fig.[4[b]), do 
not show a response any more between 3072 < t < 4096. 
The increased complexity between yc and its noisy coun- 
terpart yo is also apparent; scales 7 to 13 show a more 
irregular pattern between 4096 < t < 6144 than between 
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1024 ^ t < 3072. Finally, for the white-noise regime ys 
between 6144 ^ t < 8192, the DWaT response decreases. 
All the modes j = 1, 13 are on average less excited than 
e.g., during regime yo, and the DWaT intensity is spread 
over these scales in an almost uniform manner. This shows 
that the DWaT correctly identifies this regime as noise. We 
remark that border effects may explain the feature seen at 
scale j = 6 between < t < 1536, namely, that the magni- 
tude of that mode increases and decreases until it stabilizes 
at the end of that interval. Therefore border effects may 
here persist for up to 1536/8192 ~ 18% of the duration of 
the signal. This is still less than the conservative estimate of 
25% we have given in ^2.2.2[ Figure|4];c) shows the DWaTIM 



Hyjj . Apart from the interval < t < 1536 where the border 
effects spoil the resullj^ the DWaTIM provides a consistent 
overall diagnostic of time-dependent complexity. Each tran- 
sition between the regimes j/a and ys, through regimes yc, 
yB and yn, is detected by the DWaTIM. For instance, the 
transition a,tt = 3072 from yc to ys corresponds to a drop of 
about 0.05 in average Hyj^ magnitude around that instant. 
We recover a DWaTIM value for regime ys of Hy^j ~ 0.19 
between 3072 < i < 4096 i.e., of about 3 • lO"'^ units larger 
than in toy model 1. This 19% difference is due to the re- 
sponse of the three low- frequency scales j = 1, 2 and 3 on 
Fig.[4]at these times; in toy model 1 these modes do not show 
any response (see Fig.[3[b] between 4096 < t < 8192). The 
performance of the DWaTIM to provide an absolute mea- 
sure of complexity for a given dynamical regime depends 
on the time series subjected to analysis and on the num- 
ber of scales included in the DWaTIM computation. In §3.2| 
we investigate further this effect for the case of TV-body or- 
bits. Finally, we illustrate the behavior of the DWaTIM in 
the white-noise limit of j/b- As mentioned earlier, the total 
wavelet energy between 6144 ^t< 8192 (see grey-shade in- 
tensities on Fig.|4|b]) is reduced with respect to the interval 
4096 < t < 6144, pointing at the difficuffy of the DWaT 
to single out privileged frequencies with a high probability 
during 6144 < t < 8192. In this limit, the DWaTIM asymp- 
totically reaches its maximum value of 1 (for the case of a 
perfect white-noise signal and an infinite-length time series; 



3.1 Binary motion: N=2 

We study independently the motion of two unperturbed, 
equal-mass binaries with respective eccentricities ei — 0.05 
and 62 — 0.95. Results for the w^j-component analysis of 
body #1 of each binary are presented in the left-hand and 
right-hand panels of Fig.[5j respectively. Figure [5|^a) shows 
the velocity time series sampled at 8192 regular time 
intervals. The binaries are integrated over 16 A'^-body time 
units, so the sampling interval A is 16/8192 = 2~®. The 
period of both binaries is 0.5, allowing a sampling rate of 
256 data points per revolution. We note that a consistent 
choice of A is of major importance to the method. If an 
orbit contains frequency components that exceed /cr (see 



f 2.1.1 1 the DWaTIM will be aliased. A reasonable selection 



of A must thus ensure that the complete dynamics of the 
system is reproduced for the time-scales of interest. 

The complexity diagnostics for the two binaries are 
given in Fig.jtjb) and (c). The DWaTIM H^-^ is extracted 
from the DWaT (Fig.|5[b]) and is shown on Fig.[5];c). The 
complexity measure oscillates in time with the same fre- 
quency as the binary (in the same way than for the yi toy 
model of Fig.|3] for t > 4096). The amplitude of these os- 
cillations is greater for the more eccentric binary. This is 
so because the larger the eccentricity of a binary the larger 
the variations in velocity. In the limit where the semi-major 
axis is sufficiently large, the velocity of both bodies at apoc- 
enter is close to zero; liltewise, the complexity of any static 
configuration would asymptotically approach zero. At apoc- 
enter we find, for the e = 0.05 binary, a DWaTIM of 
whereas, for the e — 0.95 case, we obtain a value of ~ 0.1. 
In contrast, at pericenter the velocities change rapidly, re- 
quiring a broader frequency spectrum and thus implying a 
higher complexity. Here the DWaTIM indicator increases to 
~ 0.03 in the e = 0.05 case and to « 0.5 for the e — 0.95 
binary. We note that a circular binary (e = 0) has a con- 
stant velocity modulus: for this case the complexity is zero 
throughout (for all the regions not affected by border ef- 
fects, see [2.2.2). We also remarli that the time-averaged 



see f 2.2.1 1 



DWaTIM would give approximately a constant value. This 
gives a way to quantify an external perturbation acting on 



a stable Keplerian orbit (see [ 3.2 1 



We also computed the average complexity of the Vx 



3 FEW-BODY ENCOUNTERS 

To validate the technique for the case of orbital dynamics, 
we present three applications to small N problems (A = 2, 3 
and 5). AH individual orbits were computed using the star- 



lab software environment (Portegies Zwart et al. 2001 1. Inte- 
grations are performed using individual time-steps ( |Aarseth| 
1985 1 and a fourth-order Hermite predictor-corrector scheme 



(Makino & Aarseth 19921. Time series are constructed by 



Hermite interpolation at evenly spaced time intervals and by 
projecting the orbit of each star on the three orthogonal axes 
X, y and z in the center of mass coordinate system. Stan- 
dard A-body units are used throughout ( [Heggie fc Mathie"u| 
T986|. 



* The measure should give between < t < 1024. However, it 
is still true that the signal is quantified as less complex than 
yg and that the transition from j/a to yg is picked up easily. 



binaries. For each binary, the eccentricity and the period- 
icity was randomly chosen in the intervals < e < 1 and 
2~* < Hj < 2, respectively. The corresponding results for 
the DWaTIM are shown by the horizontal lines (in blue) on 
Fig.[5|c). The solid central line displays the mean DWaTIM; 
the interval delimited by the upper and the lower dashed line 
indicates the standard deviation. The average DWaTIM is 
~ 0.22. This value can be seen as indicative for the complex- 
ity of unperturbed binaries with a periodicity Hj comprised 
within the DWaT bandwidth. 



3.2 The Pythagorean problem: N=3 



The well-studied Pythagorean configuration ( |Burrau||1913[ ) 
is a classic example of long-term complex behavior (see e.g., 
Aarseth et al.|[r994 l. The initial conditions consist of three 
particles at rest, placed at the vertices of a Pythagorean tri- 
angle. The initial conditions are given in Table [T] and the 
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Figure 5. Complexity measures of two unperturbed binaries with eccentricity e = 0.05 (left-hand panels) and e = 0.95 (right-hand 
panels), respectively. For each binary, we show the a;-component result of body #1. From top to bottom: (a) velocity time series, (b) 
discrete wavelet transform (DWaT), (d) discrete wavelet transform information measure (DWaTIM) H^i- The solid horizontal lines 
in (c) show the average over the x-, y-, and 2-components for 100 binaries with e and Uj randomly chosen between < e < 1 and 
< Hj < 2, respectively. Dashed lines are one standard deviation. 



configuration is depicted in Fig. [6] 

Orbital integrations were performed over 1500 time 
units by repeatedly re-running the initial configuration with 
a reduced per-step integration error until convergence of the 
result was reached. In this way, the final trajectories showed 
a total absolute energy error |-Eflnai — -Binitiail of less than 
■|^Q-io rpj^g analysis is restrained to the first 1024 time units 
of integration, during which we found by visual inspection 
that the system spent a comparable amount of time in triv- 
ial and in more complicated states. Once more, the time 
series comprise Q = 8192 data points. Results of the respec- 
tive time series analysis for particle ^1 and for the entire 
Pythagorean problem are shown in the left-hand and right- 
hand panels of Fig. [7] respectively. 

The Hy-coordinate of body #1 is shown in the left- 
hand Fig.[7|a). In the right-hand panel of Fig.[7|a) we illus- 




m2 



Figure 6. Pythagorean problem: initial configuration. 



© 0000 RAS, MNRAS 000, 000-000 





O 500 lOOO O 5O0 1O0O 



t t 

Figure 7. Complexity measures of the x component of particle #1 (left-hand panels) and averaged overall complexity (right-hand panels) 
of the Pythagorean problem (see Tablejl] Fig.[6]and text for further detail). From top to bottom: (a) velocity time series (left-hand panel) 
and X position time series of the 3 particles (right-hand panel), (b) discrete wavelet transform (DWaT) (left-hand panel) and average 
DWaT (right-hand panel), (c) discrete wavelet transform information measure (DWaTIM) Hy^ and average DWaTIM Ttot- The dotted 
lines on (c) are obtained by computing the measures using the complete bandwidth of the DWaT. The solid lines show the same result 
computed by excluding the contribution of scales 1, ...,4. 



Body 


mass 


X 


y z 




Vy 




#1 


3 


1 


3 











#2 


4 


-2 


-1 











#3 


5 


1 


-1 












Table 1. Pythagorean problem: initial conditions. 



trate the dynamics of the whole Pythagorean problem by 
showing the position x of all the 3 bodies. Here we see that 
the system is characterized by a highly complicated inter- 
play of the three particles, including a sequence of intermit- 
tent binary formation and disruption. In the time interval 
300 < t < 540 for instance, the single body #1 strides away 
to large distance on a smooth trajectory (upper dashed line 
in the right-hand Fig.jTfa] during that interval). Due to the 



recoil, a binary, formed of bodies #2 and #3, leaves in the 
opposite direction. The motion of the binary stars about 
their common centre of mass gives the thick line seen in the 
lower part of the right-hand Fig.[7|a) between 300 < t < 540. 

The DWaT scalogram of body #1 and the average 
DWaT of the Pythagorean problem (the average magni- 
tude of base frequency l/Hj at integration time i), com- 
puted by averaging over the 3 particles and over the x- and 
y-components, are shown in the left-hand and right-hand 
panels of Fig.jtjb), respectively. The DWaTIM for body #1, 
Hy^ , and the overall DWaTIM of the Pythagorean problem, 
Ttot (see ^2.2. are shown in Fig.[7]^c). The dotted line de- 
picts the DWaTIM as computed by taking into account all 
the scales j — 1, 13. The solid line highlights the dramatic 
improvement of the diagnostic, obtained when ignoring the 
first 4 scales {j = 1, 4) in the computation. As mentioned 
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earlier (see [2.2.2), both complexity measures cease to yield 



reliable results in the limit where the orbit is resolved over 
less than 8 complete oscillations (i.e., dynamical times). This 
motivated our decision to exclude the low-frequency DWaT 
scales in a general manner. We take the conservative choice 
of not including the j = 4 scale (i.e. the periodicities 114 
corresponding to 1/8 of the full signal). For the remainder 
of this work, we only compute and discuss this improved 
version of the DWaTIM. 

The performance of the DWaTIM can be studied with 
respect to two criteria: 1) time-localization, i.e. the ability 
of the measure to identify accurately the instants t where 
qualitative changes in the orbital dynamics occur, and 2) 
continuous quantification of complexity, i.e. the ability of 
the measure to provide a consistent evaluation of complex- 
ity in an uninterrupted manner. In order to discuss 1), we 
examine the time series of body #1 of the left-hand Fig.[7]^a). 
Let us consider the transitions from a regime where body #1 
approaches bodies #2 and #3 on a comparatively smooth 
trajectory to the more complex regime where it undergoes 
multiple close encounters with the latter two bodies. By 
visual inspection of the left-hand and right-hand panels 
on Fig.[7|a), we identify such transitions to occur at e.g., 
t ~ 200, t ~ 530 and t ~ 800. In each of these cases, we 
can see in Fig.[7|c) that the measure Hy-^ shows the corre- 
sponding transition. For example, the transition at t ~ 530 
is identified by the DWaTIM by an increase of about 0.7. 
The instantaneous close encounter between the 3 particles 
at i « 800 is also singled out. Here we find a local peak of 



0.8. Likewise, the instant t « 875 at which body #1 



forms a binary system with one of the remaining bodies is 
also recovered; the complexity shows an instantaneous peak 
at that instant and then gradually decreases from Hy^ ~ 0.6 
at t « 875 to a mean value of Hy^ « 0.25 at t = 1024. The 
accuracy with which the epoch of a transition is resolved 
depends on the frequency bands involved in the transition. 
As mentioned earlier, the DWaT resolution in time increases 
by a factor of 2 when shifting from scale j to j + 1. High- 
frequency transitions can therefore be localized more accu- 
rately in time than low-frequency transitions. The maximum 
latency of the DWaTIM is obtained by studying the partic- 
ular case of a transition involving the two lower scales j = 5 
and 6. The maximum latency then corresponds to the reso- 
lution at scale j — 6, i.e. 1024/2^ = 32 time units. This sets 
an upper limit to the error on the DWaTIM performance in 
time-localization. In conclusion, qualitative changes of or- 
bital dynamics can be accurately localized in time by the 
DWaTIM indicator. 

Concerning 2) the continuous quantification of complex- 
ity, we argue that, broadly speaking, the Hy^ of the left-hand 
Fig.[7|c) reproduces the complexity of the frequency spec- 
trum obtained by the DWaT in Fig.jTJb) in a consistent 
manner. The relative complexity of the time series is quan- 
tified continuously in time with respect to the two limiting 
cases of a single base-frequency sinusoid with Hy^ — and a 
white noise signal of Hy^ — 1. Excluding the low-frequency 
scales j = 1,...,4 has also some bearings on the reliability 
of the diagnostics. For instance, when these frequencies are 
included, the DWaTIM produces inconsistent results in the 
interval 300 ^ t < 540 (see dash-dotted line in Fig.[7|c]). 
The complexity in the case of unperturbed motion of body 
#1 is then estimated to be higher than for the case where 



V3 = -V2 



T 




M 



Figure 8. Caledonian 5-body problem: initial configuration. 

the body is bound in a binary (compare to the dash-dotted 
line in the interval 875 < t < 1024). 

We also investigated the effects of analyzing the mo- 
tion of body #1 through a different time series on the Hy^ 
complexity diagnostic. This additional consistency check al- 
lows to measure DWaT border effects for the case of the 
Pythagorean problem and to examine the extent to which 
the DWaTIM is able to provide an absolute measure of com- 
plexity. To compare two time series we constructed a new 
one by taking the first half of the Q = 8192 points signal 
of Fig.lTFa), i.e. up to t = 512. This gave us a new time 
series Vy^{tq) of Q' = 4096 data points. The DWaT anal- 
ysis was then performed over a different frequency domain 
than for the original Q = 8192 {t = 1024) case. The time 
series v'y-^{tq) had also a different end value than its parent 
time series Vy^{tq), i.e. v'y^{Uo96) ^ Vyi{ts,W2). The ampli- 
tude of DWaT border effects were therefore different than 
for the Q = 8192 case (see |2.1.2[). We computed the DWa- 



TIM H'y-^ and compared it with Hy-^ shown on Fig.[7|c) in 
the interval Q <t < 512. Except for values right at the edge 
of Fig.[7|c), the results were nearly identical. In particular, 
the rms residual between the two measures was ~ 0.095 be- 
tween < t < 128 (corresponding to the first 25% of the 
signal DyJ and < 10"^ between 128 < f < 512. 

We obtain similar results for the time series of the other 
two bodies #2 and #3, and for the a;-component. For that 
reason we omitted to display those results. The overall DWa- 
TIM Ttot of the Pythagorean problem (shown on Fig.[7[c], 
right-hand panel) gives a time-resolved insight on the global 
orbital evolution obtained by integrating the initial condi- 
tions of Table[l] Once more, the features seen in the curve of 
Ttot match those in real space of Fig.[7]^a), right-hand side 
panel. We further discuss the utility of this general measure 
in the next section 93.31 



3.3 Perturbed N=5 Caledonian configuration 

Next we perturb a configuration of the planar Caledonian 
symmetric four-body problem (hereafter CP; |Szell et al.| 
2004). Initial conditions are shown in Fig.js] and in Table 
2] The CP is set up as described in Figure 1 of ISzell et| 
ah] ( |2004[ ). We specialize to the case where all bodies have 
masses equal to 1. For this configuration, the unperturbed 
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Body 



#1 


1 


0.55 





0.27 





0.23 


0.00 


#2 


1 


0.30 





0.27 





-2.42 


0.00 


#3 


1 


-0.30 





0.27 





2.42 


0.00 


#4 


1 


-0.55 





0.27 





-0.23 


0.00 


#5 


1.5 


0.00 





-0.72 





0.00 


0.03 



Table 2. Caledonian 5-body problem: initial conditions. 
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Figure 9. Overall eomplexity of the Caledonian problem (see 
also Table [2] and Fig.jsJ. From top to bottom: (a) position time- 
series of the 5 particles, (b) averaged discrete wavelet transform 
(DWaT), (c) overall information measure Ttot. 



tional interplay between the five particles is difficult to follow 
by visual inspection: tlie individual evolution of the 5 orbits 
will not be discussed in detail. We focus on the evolution of 
the system as a whole. The dynamics up to t ~ 43 may be 
described by roughly two qualitatively different states of mo- 
tion. The first one is sometimes referred to as "hierarchical 
interplay" (see e.g., Gemmeke et al. 2006). It is character- 
ized by a central binary and three particles orbiting at large 
radii. A clear example of this can be seen on Fig. [9] in the 
time interval 23 < t < 36. (The hard central binary forms 
at t ~ 23 as indicated by the vertical arrow). The second 
state of motion is what we call "democratic interplay" . This 
regime is characterized by multiple close encounters between 
the 5 bodies, each body contributing an approximately equal 
amount to the overall complexity of the system. Examples of 
democratic interplay are the time intervals 5 < t < 10 and 
14 < t < 23. At t « 43, the dynamics of the perturbed CP 
changes dramatically. After a close encounter around that 
time, all the particles stride away to larger distance. One 
particle immediately escapes the system on a nearly recti- 
linear trajectory (see upper dashed line on Fig.[9[a] at that 
instant). The remaining bodies stay close to each other until 
t ~ 46, when they are subjected to a further close encounter 
and another body is ejected (see lower solid line). The sub- 
sequent motion is relatively smooth. One particle remains 
(dashed line) and orbits around a hard binary (indicated by 
the horizontal arrow). 

The overall complexity of the perturbed CP is shown by 
the DWaTIM in Fig.|9][c). Let us consider two time intervals, 
prior to and after t = 43. When t < 43 the DWaTIM Ttot 
is roughly constant and it is interesting to observe that the 
local minima seen in that time interval can be found during 
the more quiescent hierarchical regimes, such as for exam- 
ple at i « 32. The two peaks observed at t « 43 and « 46 
arise from the two independent high-energy encounters that 
we have described earlier taking place at those times. When 
the final binary forms and the system dissolves the motion 
becomes unmistakably less complex and consequently Ttot 
has a lower value of ~ 0.2 on average for all time t > 48. 



CP consists of two stable binaries evolving in the x — y plane 
around their common center of mass. We now perturb the 
CP configuration with a fifth body of mass = 1.5 (body #5 
in Table|2|. The perturber approaches the x — y plane as it 
moves towards the center of mass of the two binaries on a 
time scale of ~ Stt/5. At that moment the distance between 
body #5 and the center of mass is ~ 0.15 i.e., about the 
separation of the two binaries. Therefore the CP is strongly 
perturbed from t ~ 5 onwards. The subsequent evolution is 
characterized by repeated interactions including the forma- 
tion of hierarchical systems with several particles orbiting 
around a hard central binary. 

The outcome of the complexity analysis is summarized 
in Fig. [9] Figure |9| a) gives the x-position time series of the 
five particles; panels (b) and (c) show the average DWaT 
and Ttot (similarly to the right-hand Figs.jTJa], [b] and [c]). 
Figure [9|a) gives an overview of the motion for the entire 
integration time of 64 time units. (Note that the motion of 
particle #5, plotted with a dotted line on Fig.[9[a], is difficult 
to disentangle and is hardly visible.) The intricate gravita- 



4 EQUAL MASS N=256 PLUMMER SPHERE 

We apply the DWaTIM technique to a self-gravitating spher- 
ical polytrope of index n = 5 (Plummer 1911 Binney & 



|Tremaiiie||1987[ §4.4.3). The phase-space distribution func- 
tion of this system is a power-law of e, 



24^2 

77r3 GSM* 



(26) 



where e = ^v'^ + ${r) < is the mechanical energy per 
unit mass, v the three-dimensional velocity, $ the potential 
which is a function of the radius r only. Note that = 
by construction whenever e > so that only mass elements 
bound by gravity are considered. Given a value of the gravi- 
tational constant G, the two free parameters Af and R define 
the total system mass and a reference unit of length, respec- 
tively. Integrating over all velocities at constant radius 
yields the mass density g at that radius. 



g{r) = 



V-2*('-) 



'iTTj-{e)v dv = 



_3_ M 



1 + 



R 



-5/2 



(27) 
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Figure 10. From bottom to top: core radius rc, 20%, 30%, 40%, 
50%, 60%, 70%, 80% and 90% Lagrangian radii of the N = 256 
equal-mass Plummer sphere. The soUd vertical line at t = 113 
marks the formation of the first hard binary in the system. The 
dashed vertical lines indicate the onset and the end of a notable 
core expansion between 165 ^ i ^ 180. 



where we have used d^i; = Anv^dv valid for an isotropic 
velocity field and where we have substituted for <l?(r) by 
solving Poisson's equation. Note that the length R defines 
the radius of a uniform-density core {g[r << R] « con- 
stant). An A'^ = 256 particle representation of Eq. ( |27[ ) is ob- 
tained by random-sampling the mass density to assign three- 
dimensional positions. All bodies have mass rrii = M/N , 
where i = 1, ...,N. The particles' energy is attributed simi- 
larly from Eq. (|26|) from which we compute the squ are veloc- 

= 2(£ 



Aarseth, 



ity = 2(£ - $[r ]) as in the standard method of 
[Henon fc Wielen] ( |1974| |. In those circumstances the total 
kinetic and gravitational energies satisfy the virial theorem 
of equilibrium systems. 

The nominal dynamical time is = 2R/a where is 
the mean squared velocity averaged over mass up to the half- 
mass radius. We found = 0.392...GM/R ~ nGM/8R. 
The two-body relaxation time is conveniently defined as 



t, = Ntd/\riOAN (see |Binney fc Tremaine||1987[ ). If we set 
Heggie-Mathieu computational units with G = M = 1 and 
R = 1.1781 such that the total binding energy E — —1/4, 
we obtain td ~ 3.46 so that for A'' = 256 bodies the relax- 
ation time ti « 191 time units. For low-A'^ systems such as 
this one, the diffusion of kinetic energy leads to core-collapse 
on roughly that time scale, when q becomes singular at the 
center. It is around that time that hard binaries form and 
energy exchanges between bodies is at its most extreme. We 
therefore evolved the system for a total time of 256 time 
units to ensure that the cluster reaches core-collapse and 
contrast pre- and post-collapse evolution. The equations of 
motion were integrated with no softening of the potential. 



The velocities of all bodies are stored at regular time 
intervals and the DWaTIM is computed in a post-simulation 
analysis. The sampling interval. A, should be sufficiently 
small to capture the dynamics in the dense core fully. A naive 
strategy would be to store all data at all integration time 
steps. For the configuration adopted here we have observed 
that the time steps may become as small as 2~^^ ~ 2.32 • 
10~^°; unfortunately a sampling rate of that order would 
translate to a prohibitively high volume of data storage even 
for modest values of A^. Instead, we note that the dynamical 
time td in principle sets a standard from which to pick an 
adequate value of A. But because the central region becomes 
ever denser during evolution, td ~ 3.46 set from the initial 
configuration would give no useful reference. 

Instead, we explore the evolution in time of the mass 
profile of the system and look for a minimum core length and 
velocity dispersion inside that core. The density-averaged 
core radius defines a quantity which monitors the rise of 
the central density (see e.g. von Hoerner 1960, Casertano 



& Hut 1985). From Eq. (271 we find for the initial configu- 
ration rc — R/2 ~ 0.59 (enclosing 13% of the mass, or 33 
bodies). Figure [To| displays the time-evolution of rc and sev- 
eral Lagrange radii. The core radius rc always appears at the 
bottom on Fig. |10[ Up to f ~ 113 units, rc decreases on the 
mean. At that time, the first hard binary of binding energy 
E < -lOOkT formf[f] That event is marked with a vertical 
full line on Fig.[To] Note that rc along with the 20% and 
30% Lagrange radii increase significantly from t ~ 165 on- 
wards. This phase of rapid expansion indicates the on-set of 
post-collapse evolution for that simulation. The minimum 
value of rc — 0.05 occurs at f ~ 165 units and encloses 
3% of the total mass (8 bodies). We compute a dynamical 
time td,c — 0.05 for the core at that time, a factor ~ 70 
smaller that the value computed from the initial conditions. 
An orbit confined to the core is adequately sampled with 
five points or more and hence we set A = 0.01 for the com- 
plexity analysis. The time-resolution of any features seen in 
the diagnostic of complexity is therefore 2A — 0.02 units, 
and any binary formed through dynamical evolution is well 
sampled provided its binding energy E > —2.75kT. 



4.1 Individual orbits 

To illustrate the differences in orbital complexity a star may 
show between the initial time and the moment of core expan- 
sion, we graph in Fig.[TT]the x — y projections of 2 individual 
orbits during two windows of four time units, the first run- 
ning from t — to 4: (top panels), the second running in the 
interval t — 160 to 164 (bottom panels). We picked two stars 
that happened to orbit within the core in each time interval 
(the circles on the figure indicate rc at the times shown). 
In both the cases displayed the orbit starts off smooth and 
regular, but traces a much more intricate pattern later on. 

These trends can be identified in a graph of the DWa- 
TIM for these and four other orbits as displayed on Fig. |12| 
The figure shows vs time in the main frames and the 
DWaTIM as the top inset frame for each case. Note the 



The energy scale kT is defined by the condition that the to- 
tal stellar kinetic energy of the system, excluding internal binary 
motion, is 3/2NkT. 
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Figure 11. Two stars evolving in tlie core of tiie A'^ = 256 equal-mass Plummer sphere: (a) the x — y projections of particles #6 and 
#182, evolving in the core over 4 time units starting at t = 0, (b) their x — y projections starting at t = 160. The outer dashed circle 
in (a) shows the core radius Vc of the cluster at 4 = 0. The inner circle in (a) depicts the core radius at t = 160 as shown in (b) and 
illustrates the difference in scale between the two figures. 



change of scales: the DWaTIM is plotted against the scale 
shown at the right-hand side of each figure. Focusing on the 
top two panels on Fig. |12[ we can identify the more com- 
plex phases around 160 < t < 164 of the orbits displayed on 
Fig. 11 b) as local peaks in the DWaTIM during this time 



post-collapse core expansion, respectively (cf. also Fig. 10 1 



interval. It is clear that stars set on regular orbits initially 
can show more complex behavior at later times. The oppo- 
site is also possible, as we illustrate in Fig. [T2]with four more 
orbits also orbiting in and out of the core region. These and 
many others not shown here are typical of the wide variety 
of DWaTIM spectra: some orbits show rapid fluctuations in 
Vx and yield a rather broad band of base frequencies (es- 
pecially particles #6 and #182). Other trajectories have a 
narrower spectrum of frequencies (e.g., particle #128). 



4.2 Global behavior of the Plummer sphere 

Fig.[T3] shows as a function of time the cumulative num- 
ber of stars whose DWaTIM Ti (i = 1,...,A'') exceeds a 
given threshold. The uppermost horizontal line on the fig- 
ure denotes the total number of stars. The broken curves 
are for (from top to bottom) Ti > 0.2, 0.4, 0.5, 0.6 and 0.7. 
Once more, the vertical solid and dashed lines indicate the 
time when a hard binary first formed and the interval of 



Broadly speaking, the orbital complexity decreases on the 
mean with time. The uppermost lines showing > 0.2 and 
Ti > 0.4 decrease monotonically save for small localised 
fiuctuations. In all the cases displayed, a pronounced drop 
in complexity is seen throughout the phase of core expan- 
sion, at times 160 ^ t < 185. For instance, at the beginning 
of the expansion, about 110 stars have Ti > 0.4 whereas 
at t m 180 only 75 stars reach that level. The situation is 
similar for higher-threshold curves. After the formation of a 
hard binary but prior to post-collapse expansion, i.e. in the 
interval 113 ^ t < 165, the complexity levels off or increases 
slightly with time. For example, the number of stars with 
Ti > 0.5 increases from about 40 to approximately 55 stars 
in that interval. Hence, orbits that are already relatively 
complex at the time of binary formation yield a DWaTIM 
of even larger amplitude up to core-collapse and the on-set 
of the expansion phase. This trend is also found in the curve 
Ti > 0.6 and to a lesser extent in the Ti > 0.7 curve. It may 
be important to note that the global results of Fig.[T3] are 
not in contradiction with the apparent trend of increasing 
complexity depicted by the two stars of Fig. |ll| Stars #6 and 
#182 of that figure are both part of the 8 stars that happen 
to reside in the core around t ^ 165. At these times the stars 
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Figure 12. Time evolution of 6 individual orbits in an N = 256 equal-mass Plummcr sphere model. Each figure shows the complexity 
diagnostic obtained with respect to the x-coordinate of particle i, and comprises the velocity time series Vx^ (solid line in red, plotted 
against the left-hand axis) and the DWaTIM Hx^ (solid blue line in the upper inset, plotted against the 0—1 scale shown on the 
right-hand axis). 
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Figure 13. Cumulative number of stars of the A' = 256 equal- 
mass Plummer sphere with a DWaTIM {i = l,...,N) larger 
than (from top to bottom) 0.2, 0.4, 0.5, 0.6 and 0.7, respectively. 
The uppermost horizontal line shows the total number of particles 
A'^ = 256. The solid vertical line at t = 113 marks the formation 
of the first hard binary in the system. The dashed vertical lines 
indicate the phase of core expansion (165 ^ 180). 



have a high DWaTIM with values of Ts and T182 > 0.6 (see 
also the two upper panels of Fig. |12[ ). Stars #6 and #182 
therefore contribute to the local peak of the Ti > 0.6 curve 
observed on Fig. [13] around t « 160. The snapshots given 
in Fig. |ll[ b) do not reflect a progressive increase in orbital 
complexity of the Plummer sphere between < t < 164. 
The enhanced two-body scatter attributable to the newly 
formed binary can be directly measured by the DWaT. In 
conclusion, both that event and the on-set of post-collapse 
expansion can be singled out on Fig.[T3]as a local minimum 
and a local maximum, respectively, in runs of the cumulative 
Ti (as examplified, for instance, by the curve of Ti > 0.5). 

A scalogram of the DWaT averaged over all = 256 
Plummer sphere particles is shown on Fig. |14| The dark 
shade illustrates base frequencies of high amplitude. 



whereas Fig 
j 



white means zero amplitude (see also ^3.2 and S3.3I. Figure 
T4|a) brushes a global picture for all scales j = 1, 15, 
b) depicts only the high-frequency scales 
11, 12, and 13. Some modes are growing in intensity 
and are then fading away after t ~ 165 units when the core 
starts to expand. A particularly good example is the scale 
j — 6 which reaches progressively higher amplitude in the in- 
terval < t < 160 before fading away during the expansion 
phase of the inner volume of the sphere. The onset of ex- 
pansion triggers high-frequency modes (cf. Fig. 14 b]) which 



reflect the evolution toward more anisotropic radial orbits. 
Radial anisotropy implies more eccentric motion relative to 
the centre of mass and an enhanced spectrum of frequencies, 
cf. Fig.[5l 

An illustration of the global measure of complexity 





Figure 14. Average discrete wavelet transform (DWaT) of the 
N = 256 Plummer sphere as a function of time. The averaging 
is performed for a total of 768 time series, i.e. over all particles 
and over the three velocity components. The figure illustrates 
the average amplitude of base frequency l/Tlj: (a) the complete 
DWaT scalogram (scales j = 1,...,15); (b) zoom on the high- 
frequency scales j = 11, 12, and 13. 



15 



The solid and dashed 



Ttot — T^i/N is shown on Fig. 
vertical lines denote the same transitional phenomena as in 
all preceding figures. On the whole, the curve of the DWa- 
TIM depicts the same global decrease in complexity as ob- 
served on Fig. |13[ A close inspection of Fig.[T5]suggests three 
different regimes in the evolution of the DWaTIM indicator 
Ttot '■ 

- Regime 1 takes place between < i < 165. In the 
interval < t < 20, we observe a rapid drop in Ttot which 
has no clear origin in a physical phenomenon (binary for- 
mation, core-collapse, etc). The equilibrium cluster has a 
dynamical time — 3.5 and hence the elapsed time t = 20 
is « 8td or eight orbital times, a further illustration of bor- 
der effects arising from a too-short time of integration. For 
t > 20, the curve decreases on average until t ~ 110 when 
a slower but constant decline sets in. The DWaTIM Ttot 
decreases from « 0.38 at t ~ 20 to Ttot « 0.35 at t ~ 165; 
the formation of soft- and hard-binaries therefore has little 
impact on the global complexity of the system. That being 
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said, it is difficult to disentangle the apparent trend of a 
drop in complexity during regime 1 because, first of all, the 
rapid drop early on is attributable to border effects; and 
second, the time scale over which the trend becomes sig- 
nificant is comparable to the two-body relaxation time for 
that system. The data on Fig. [15] may be best understood 
in the light of Fig. |13[ In that figure stars of low complexity 
show a decrease of between < t < 110 (see the two 
upper curves Ti > 0.2 and > 0.4) whereas the number of 
stars on orbits giving a diagnostics of high-complexity re- 
mains approximately constant. This argues against border 
effects reaching beyond t ~ 20. The constant evolution to- 
wards more orbits of low-complexity diagnostics drive the 
trend of Ttot observed on Fig. |15| a statement that the sys- 
tem suffers from coUisional effects when stars are shifted to 
higher-energy long-period orbits by two-body collisions. In- 
spection of the rapid evolution of the core-radius comforts 
this view. 

- Regime 2 begins around t ~ 165 when the system en- 
ters post-collapse and the central volume expands system- 
atically. This triggers a sudden drop in the DWaTIM from 
« 0.35 to ^ 0.27 at t ~ 180. 

- Regime 3 begins (loosely speaking) at t ~ 180 units 
when the run of Ttot resumes a slow decrease on the average. 
This contrast in DWaTIM orbital complexity between the 
pre- and the post-collapse evolution has been seen in a set 
of ten N = 256 test simulations we performed with differ- 
ent random seeds. An attempt to obtain ensemble averaged 
results for this set of simulations proved fruitless owing to 
large scatter e.g., in the core-collapse time between two in- 
dividual runs. 




Figure 15. Overall DWaTIM complexity Ttot for an = 256 
equal-mass Plummer sphere. The curve is the average of over 
all N particles and over the three velocity time series in each case 
(the X—, y— and z— components). The solid vertical line at t = 
113 together with the right-most arrow mark the formation of the 
first hard (stable) E/kT < —100 binary in the system. The left- 
most arrow marks the time when the first soft (unstable) E/kT < 
— 10 binary formed. The dashed vertical lines between 165 ^ 
180 bracket the period of expansion of the inner volume. 



5 DISCUSSION 



We presented a method to compute the time-dependent or- 
bital complexity in A'^-body simulations. The gravitational 
A-body problem is described by the 3N second order ordi- 
nary differential equations of Eq. ([l]). We extract a discrete 
wavelet transform information measure (DWaTIM) from the 
velocity time series of the individual particles of the A-body 
simulation. We apply the technique to several few-body 
problems and to a larger N = 256 particle simulation. The 
method captures the time-dependent changes in the dynam- 
ics of three- and five-body systems and furthermore quan- 
tifies orbital complexity continuously in time. For example, 
we recovered and quantified the dynamically more complex 



phases of the well-studied Pythagorean problem (see S3.2I 
as well as the complex dynamics of a perturbed Caledonian 



configuration (see f3.3l. We also applied the method to a 
set oi N = 256 equal-mass Plummer spheres (see We 
found that, on a global scale, orbital complexity decreases 
during the evolution. The occurrence of core-collapse and 
the subsequent core expansion causes a considerable drop 
in overall DWaTIM complexity Ttot- Furthermore, we ob- 
served that the complexity of individual orbits with a DWa- 
TIM Ti > 0.5 (i = 1,...,N) at the instant of core-collapse 
tends to increase until the occurrence of core expansion. 

We opted for a DWaTIM implementation that allowed 
1) to identify qualitative changes in orbital dynamics in a 
quasi-instantaneous manner and 2) to provide a continu- 
ous quantification of time-dependent complexity on a well- 



defined gauge between and 1. One should however bear 
in mind that an absolute measure of complexity of a dy- 
namical system cannot be obtained by the method. Com- 
plexity ultimately depends on the range of scales over which 
the system is studied. The DWaTIM is a band-limited mea- 
sure and therefore strongly depends on the resolution of the 
time series that is subjected to analysis. Experience drawn 
from several test cases shows that a repeat calculation with 
truncated frequency range is desirable to confirm the con- 
vergence of the diagnostics. 

Border effects are an unwanted artifact of the wavelet 
transform. In order to reduce the bias so introduced we 
opted to implement a discrete wavelet transform and to se- 
lect cubic spline functions as mother wavelets. This DWaT 
implementation avoids redundant wavelet coefficients. This 
helps in obtaining a proper measure of complexity for the 
two limiting cases of a unique base frequency sinusoid (DWa- 
TIM = 0) and of a white noise signal (DWaTIM = 1), and 
thus to recover a well-defined gauge of complexity. We ob- 
tain reliable diagnostics of complexity whenever the orbit is 
integrated over a minimum of 8 complete oscillations (see 

The method is computationally inexpensive. For exam- 



ple, the analysis of a time series of length Q 



data 



points takes approximately 1 second on a Pentium IV 2.4 
GHz workstation with 1.2 GB RAM memory. It is then pos- 
sible to obtain a complexity diagnostic of an A = 10^ body 
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system in about 80 processor-hours. A parallel implementa- 
tion of the scheme is currently under development. 
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